rm(list=ls())
library(foreign)
library(stargazer)
library(lfe)
library(car)
library(ggplot2)
library(readstata13)
library(tidyr)
library(dplyr)
library(readstata13)
setwd('/Users/Jacque/Documents/Mine/Research & Study/Rochester/Projects/Dissertation/The industrial organization origin/Data_Submission_New/Final')
#########################################################################################################################################
#                                                Part I: STRIKE                                                                        #
#########################################################################################################################################
combine <- read.dta('combine.dta')
combine = subset(combine, !is.na(approve_agg))

combine$approve_agg_binary = 0
for (i in 1:length(combine$cty)){
  if (combine$approve_agg[i] == 1){
    combine$approve_agg_binary[i] = 1
  }
}

#Table C.2
t1_non = subset(combine, democracy == 0 & approve_agg_binary == 1)
t2_non = subset(combine, democracy == 0 & approve_agg_binary == 0)
a = t.test(t1_non$domestic2, t2_non$domestic2, alternative = "greater")
t1_dem = subset(combine, democracy == 1 & approve_agg_binary == 1)
t2_dem = subset(combine, democracy == 1 & approve_agg_binary == 0)
b = t.test(t1_dem$domestic2, t2_dem$domestic2, alternative = "greater")

#Figure 7
marginal1.coefficient = c(mean(t1_non$domestic2)- mean(t2_non$domestic2,na.rm=TRUE), mean(t1_dem$domestic2)- mean(t2_dem$domestic2,na.rm=TRUE))
marginal1.ci_low = c(a$conf.int[1],b$conf.int[1])
df <- data.frame(PolityIV =0:1, FDI_Restriction = marginal1.coefficient, L = marginal1.ci_low)
df$PolityIV = as.character(df$PolityIV)
df$PolityIV[1] = "Dictatorship"
df$PolityIV[2] = "Democracy"
colnames(df)[colnames(df)=="PolityIV"] <- "RegimeType"
require(ggplot2)
ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = RegimeType, ymax = FDI_Restriction, ymin = L), width=0.1, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = RegimeType, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "General Strike Difference", x = "Regime Type") +
  ggtitle("General Strikes") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))

#Table C.3
s1_non = subset(combine, polity4 < 5 & approve_agg_binary == 1)
s2_non = subset(combine, polity4 < 5 & approve_agg_binary == 0)
t.test(s1_non$domestic2, s2_non$domestic2, alternative = "greater")

s1_dem = subset(combine, polity4 >= 5 & approve_agg_binary == 1)
s2_dem = subset(combine, polity4 >= 5 & approve_agg_binary == 0)
t.test(s1_dem$domestic2, s2_dem$domestic2, alternative = "greater")

#Table 3 Columns 1 and Table C.6 Columns 1
strike1 <- felm(domestic2 ~ approve_agg + bank + currency + reschedule + signed_imf + 
                  log(gdp_pc) + trade_gdp + democracy + urban_pop + approve_agg:democracy |cty+year|0|cty, data = combine)
summary(strike1)

#Figure 8
marginal1.coefficient = c(strike1$coefficients[1], strike1$coefficients[1] + strike1$coefficients[10])
marginal1.variance = c(strike1$vcv[1,1], strike1$vcv[1,1] + strike1$vcv[10,10] + 2*strike1$vcv[1,10])
marginal1.sd = sqrt(marginal1.variance)
marginal1.ci_up = c(sqrt(marginal1.variance[1])*qnorm(0.95) + marginal1.coefficient[1],
                    sqrt(marginal1.variance[2])*qnorm(0.95) + marginal1.coefficient[2])
marginal1.ci_low = c(sqrt(marginal1.variance[1])*qnorm(0.05) + marginal1.coefficient[1],
                     sqrt(marginal1.variance[2])*qnorm(0.05) + marginal1.coefficient[2])
df <- data.frame(PolityIV =0:1, FDI_Restriction = marginal1.coefficient, L = marginal1.ci_low, U = marginal1.ci_up)

require(ggplot2)
ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = PolityIV, ymax = U, ymin = L), width=0.1, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = PolityIV, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "Marginal Effect of \nGovernment Screening", x = "Democracy") +
  ggtitle("General Strikes") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))

#Table 3 Columns 2 and Table C.6 Columns 2
strike1.1 <- felm(domestic2 ~ approve_agg + bank + currency + reschedule + signed_imf + 
                    log(gdp_pc) + trade_gdp + polity4 + urban_pop + approve_agg:polity4 |cty+year|0|cty, data = combine)
summary(strike1.1)

#Figure 8
marginal1.2.coefficient = c(strike1.1$coefficients[1] -10*strike1.1$coefficients[10], 
                            strike1.1$coefficients[1] -9*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -8* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -7* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -6* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -5* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -4* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -3* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] -2* strike1.1$coefficients[10],
                            strike1.1$coefficients[1] - strike1.1$coefficients[10],
                            strike1.1$coefficients[1], 
                            strike1.1$coefficients[1] + strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 2*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 3*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 4*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 5*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 6*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 7*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 8*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 9*strike1.1$coefficients[10],
                            strike1.1$coefficients[1] + 10*strike1.1$coefficients[10])
marginal1.2.variance = c(strike1.1$vcv[1,1] + 100 * strike1.1$vcv[10,10] - 10*2*strike1.1$vcv[1,10], 
                         strike1.1$vcv[1,1] + 81 * strike1.1$vcv[10,10] -9* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 8*8 * strike1.1$vcv[10,10] -8* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 7*7 * strike1.1$vcv[10,10] -7* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 6*6 * strike1.1$vcv[10,10] -6* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 5*5 * strike1.1$vcv[10,10] -5* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 4*4 * strike1.1$vcv[10,10] -4* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 3*3 * strike1.1$vcv[10,10] -3* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 2*2 * strike1.1$vcv[10,10] -2* 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + strike1.1$vcv[10,10] -2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1],
                         strike1.1$vcv[1,1] + strike1.1$vcv[10,10] + 2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 2*2*strike1.1$vcv[10,10] + 2*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 3*3*strike1.1$vcv[10,10] + 3*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 4*4*strike1.1$vcv[10,10] + 4*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 5*5*strike1.1$vcv[10,10] + 5*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 6*6*strike1.1$vcv[10,10] + 6*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 7*7*strike1.1$vcv[10,10] + 7*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 8*8*strike1.1$vcv[10,10] + 8*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 9*9*strike1.1$vcv[10,10] + 9*2*strike1.1$vcv[1,10],
                         strike1.1$vcv[1,1] + 10*10*strike1.1$vcv[10,10] + 10*2*strike1.1$vcv[1,10])
marginal1.2.sd = sqrt(marginal1.2.variance)
marginal1.2.ci_up = marginal1.2.sd*qnorm(0.95) + marginal1.2.coefficient
marginal1.2.ci_low = marginal1.2.sd*qnorm(0.05) + marginal1.2.coefficient
df <- data.frame(PolityIV =-10:10, FDI_Restriction = marginal1.2.coefficient, L = marginal1.2.ci_low, U = marginal1.2.ci_up)

ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = PolityIV, ymax = U, ymin = L), width=0.2, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = PolityIV, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "Marginal Effect of \nGovernment Screening", x = "Polity4") +
  ggtitle("General Strikes") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))



#########################################################################################################################################
#                                                Part II: VDEM                                                                          #
#########################################################################################################################################
#Table C.4
a = t.test(t1_non$v2exbribe, t2_non$v2exbribe, alternative = "less")
b = t.test(t1_dem$v2exbribe, t2_dem$v2exbribe, alternative = "less")

#Figure 7
marginal1.coefficient = c(mean(t1_non$v2exbribe,na.rm=TRUE)- mean(t2_non$v2exbribe,na.rm=TRUE), mean(t1_dem$v2exbribe,na.rm=TRUE)- mean(t2_dem$v2exbribe,na.rm=TRUE))
marginal1.ci_low = c(a$conf.int[2],b$conf.int[2])
marginal1.ci_high = c(a$conf.int[1],b$conf.int[1])
df <- data.frame(PolityIV =0:1, FDI_Restriction = marginal1.coefficient, U = marginal1.ci_low, L = marginal1.ci_high)
df$PolityIV = as.character(df$PolityIV)
df$PolityIV[1] = "Dictatorship"
df$PolityIV[2] = "Democracy"
colnames(df)[colnames(df)=="PolityIV"] <- "RegimeType"

ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = RegimeType, ymax = U, ymin = FDI_Restriction), width=0.1, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = RegimeType, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "Corruption Difference", x = "Regime Type") +
  ggtitle("Control of Executive Corruption") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))

#Table C.5
t.test(s1_non$v2exbribe, s2_non$v2exbribe, alternative = "less")
t.test(s1_dem$v2exbribe, s2_dem$v2exbribe, alternative = "less")

#Table 3 Columns 3 and Table C.6 Columns 3
vdem1.1 <- felm(v2exbribe ~ approve_agg + approve_agg:democracy +democracy + currency + 
                  bank + reschedule + signed_imf + log(gdp_pc) + trade_gdp + urban_pop|cty+year|0|cty, 
                data = combine)
summary(vdem1.1) 

#Figure 8
marginal1.coefficient = c(vdem1.1$coefficients[1], vdem1.1$coefficients[1] + vdem1.1$coefficients[10])
marginal1.variance = c(vdem1.1$vcv[1,1], vdem1.1$vcv[1,1] + vdem1.1$vcv[10,10] + 2*vdem1.1$vcv[1,10])
marginal1.sd = sqrt(marginal1.variance)
marginal1.ci_up = c(sqrt(marginal1.variance[1])*qnorm(0.95) + marginal1.coefficient[1],
                    sqrt(marginal1.variance[2])*qnorm(0.95) + marginal1.coefficient[2])
marginal1.ci_low = c(sqrt(marginal1.variance[1])*qnorm(0.05) + marginal1.coefficient[1],
                     sqrt(marginal1.variance[2])*qnorm(0.05) + marginal1.coefficient[2])
df <- data.frame(PolityIV =0:1, FDI_Restriction = marginal1.coefficient, L = marginal1.ci_low, U = marginal1.ci_up)

ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = PolityIV, ymax = U, ymin = L), width=0.1, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = PolityIV, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "Marginal Effect of \nGovernment Screening", x = "Democracy") +
  ggtitle("Control of Executive Corruption") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))


#Table 3 Columns 4 and Table C.6 Columns 4
vdem1.2 <- felm(v2exbribe ~ approve_agg + approve_agg:polity4 + polity4 + currency + 
                  bank + reschedule + signed_imf + log(gdp_pc) + trade_gdp+ urban_pop|cty+year|0|cty, data = combine)
summary(vdem1.2) 

#Figure 8
marginal1.2.coefficient = c(vdem1.2$coefficients[1] -10*vdem1.2$coefficients[10], 
                            vdem1.2$coefficients[1] -9*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -8* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -7* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -6* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -5* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -4* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -3* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] -2* vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] - vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1], 
                            vdem1.2$coefficients[1] + vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 2*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 3*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 4*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 5*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 6*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 7*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 8*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 9*vdem1.2$coefficients[10],
                            vdem1.2$coefficients[1] + 10*vdem1.2$coefficients[10])
marginal1.2.variance = c(vdem1.2$vcv[1,1] + 100 * vdem1.2$vcv[10,10] - 10*2*vdem1.2$vcv[1,10], 
                         vdem1.2$vcv[1,1] + 81 * vdem1.2$vcv[10,10] -9* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 8*8 * vdem1.2$vcv[10,10] -8* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 7*7 * vdem1.2$vcv[10,10] -7* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 6*6 * vdem1.2$vcv[10,10] -6* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 5*5 * vdem1.2$vcv[10,10] -5* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 4*4 * vdem1.2$vcv[10,10] -4* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 3*3 *vdem1.2$vcv[10,10] -3* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 2*2 * vdem1.2$vcv[10,10] -2* 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + vdem1.2$vcv[10,10] -2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1],
                         vdem1.2$vcv[1,1] + vdem1.2$vcv[10,10] + 2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 2*2*vdem1.2$vcv[10,10] + 2*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 3*3*vdem1.2$vcv[10,10] + 3*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 4*4*vdem1.2$vcv[10,10] + 4*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 5*5*vdem1.2$vcv[10,10] + 5*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 6*6*vdem1.2$vcv[10,10] + 6*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 7*7*vdem1.2$vcv[10,10] + 7*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 8*8*vdem1.2$vcv[10,10] + 8*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 9*9*vdem1.2$vcv[10,10] + 9*2*vdem1.2$vcv[1,10],
                         vdem1.2$vcv[1,1] + 10*10*vdem1.2$vcv[10,10] + 10*2*vdem1.2$vcv[1,10])
marginal1.2.sd = sqrt(marginal1.2.variance)
marginal1.2.ci_up = marginal1.2.sd*qnorm(0.95) + marginal1.2.coefficient
marginal1.2.ci_low = marginal1.2.sd*qnorm(0.05) + marginal1.2.coefficient
df <- data.frame(PolityIV =-10:10, FDI_Restriction = marginal1.2.coefficient, L = marginal1.2.ci_low, U = marginal1.2.ci_up)

ggplot() + 
  geom_errorbar(data = df, mapping = aes(x = PolityIV, ymax = U, ymin = L), width=0.2, size=0.5, color="blue") + 
  geom_point(data = df, mapping = aes(x = PolityIV, y = FDI_Restriction), size = 4, shape=21, fill="white") + 
  labs(y= "Marginal Effect of \nGovernment Screening", x = "Polity4") +
  ggtitle("Control of Executive Corruption") + geom_hline(yintercept=0, color = 'black') + 
  theme(plot.title = element_text(hjust = 0.5))

#########################################################################################################################################
#                                                Part III: Wage                                                                       #
#########################################################################################################################################
#Table C.6, Column 5
wage1 <- felm(log(Value) ~ approve_agg + approve_agg:polity4 + polity4 + currency + 
                bank + signed_imf + log(gdp_pc) + trade_gdp+ urban_pop|cty+year|0|cty, data = combine)
summary(wage1)

#########################################################################################################################################
#                                                Part IV: Tax                                                                          #
#########################################################################################################################################
#Table C.7
tax1.1 <- felm(tax_nr_corp ~ approve_agg + approve_agg:democracy +democracy + currency + 
                  bank + reschedule + signed_imf + log(gdp_pc) + trade_gdp + urban_pop|cty+year|0|cty, 
                data = combine)
summary(tax1.1) 

tax1.2 <- felm(tax_nr_corp ~ approve_agg + approve_agg:polity4 + polity4 + currency + 
                 bank + reschedule + signed_imf + log(gdp_pc) + trade_gdp + urban_pop|cty+year|0|cty, 
               data = combine)
summary(tax1.2) 

#########################################################################################################################################
#                                                Part IV: Korea                                                                          #
#########################################################################################################################################
#Figure 6
mean_wage = read.csv("mean_wage_Korea.csv")
plot(x = mean_wage$year, y = mean_wage$real_mean_wage,type = 'b', main = 'Monthly Wage in LCU', xlab = "Year", ylab = "Monthly Wage in  LCU", col="blue")
abline(v = 1987, col = 'red', lty=c(1,2), lwd=c(1, 3))

KOR <- read.csv('KOR.csv')
plot(x = KOR$year, y = KOR$ws, type = 'b', main = 'Labor Share', xlab = "Year", ylab = "Labor Share", col="blue")
abline(v = 1987, col = 'red', lty=c(1,2), lwd=c(1, 3))